For example, when performing an image reconstruction in an X-ray CT (Computed Tomography) apparatus as a tomography apparatus, filtered backprojection (FBP: Filtered Back Projection) method has been conventionally used as a standard image reconstruction method. On the other hand, in recent years, research and practical realization of the image reconstruction using an iterative method is progressing due to, e.g., improved performance of computers. In the case of using an iterative method, to reduce artifacts due to various factors, it is possible to reflect complicated physical models, prior knowledge, etc., and therefore various methods have been proposed so far (see Patent Document 1, Patent Document 2, Non-Patent Document 1).
Such a method can be considered as an iterative method based on objective function maximization. In this method, a reconstructed image is obtained by maximizing the objective function F represented by the following Formula (1).F(μ,y)=D(μ,y)+βR(μ)  Formula (1)
Here, μ in the aforementioned Formula (1) is a reconstructed image vector, and y is projection data. D is called “data term” and the like, which indicates the degree of conformity with measurement data, and is defined by likelihood calculated from measured projection (actually measured projection data obtained by an X-ray detector) and estimated parameters (the image estimated by the aforementioned Formula (1)). Since μ and y are vectors, they are written in bold.
In addition, R is generally called “penalty term” or the like and reflects the validity of the estimated parameter (estimated image). In this specification, hereinafter, R will be referred to as “validity term” for the sake of convenience. Further, β is a coefficient that controls the strength of the validity term R, which is determined empirically.
In the actual calculation in the aforementioned Formula (1), an optimization algorithm such as a steepest descent method and a Newton method is used. Further, in order to avoid falling into a local solution, a combinatorial optimization method such as a genetic algorithm and an annealing method may sometimes be incorporated. When a steepest descent method is used as an optimization algorithm, the update formula of the reconstructed image by the aforementioned objective function is represented by the following Formula (2).
      Formula    ⁢                              ⁢                            (    2    )                                                                                            µ                                      n                    +                    1                                                  =                                ⁢                                                      µ                    n                                    +                                      α                    ×                                          ∇                                              F                        ⁡                                                  (                                                      µ                            ,                            y                                                    )                                                                                                                                                                                            =                                ⁢                                                      µ                    n                                    +                                      α                    ×                                          ∇                                              D                        ⁡                                                  (                                                      µ                            ,                            y                                                    )                                                                                                      +                                      α                    ×                    β                    ×                                          ∇                                              R                        ⁡                                                  (                          µ                          )                                                                                                                                                                              (                                  ⁢                                                  ⁢            2                    )                    
α in Formula (2) is called a step size and has a role of controlling the magnitude of the update amount in the gradient direction. Further, ∇ in Formula (2) is a gradient, which is a partial differential with respect to the estimation parameter (reconstructed image). The update formula in the jth pixel is represented by the following Formula (3).
                              μ          j                      n            +            1                          =                              μ            j            n                    +                      α            ×                                          ∂                                                                                              ∂                                  μ                  j                                                      ⁢                          D              ⁡                              (                                  µ                  ,                  y                                )                                              +                      α            ×            β            ×                                          ∂                                                                                              ∂                                  μ                  j                                                      ⁢                          R              ⁡                              (                µ                )                                                                        Formula        ⁢                                  ⁢                  (          3          )                    
In the aforementioned Formula (3), α and β are assumed to not depend on the number of iterations and the pixel position, and R is assumed to not depend on the number of iterations. On the other hand, they may also be assumed to depend on the number of iterations and the pixel position. For example, the aforementioned Non-Patent Document 2 is directed to an image reconstruction method in which α depends on the number of iterations (indirectly also depends on the pixel value which is being reconstructed), and the aforementioned Patent Document 2 is directed to an image reconstruction method in which α and β depend on the pixel position.
In this way, when α, β, and R are assumed to depend on the number of iterations and the pixel position, the aforementioned Formula (3) may be represented by the following Formula (4).
                              μ          j                      n            +            1                          =                              μ            j            n                    +                                    α              j              n                        ×                                          ∂                                                                                              ∂                                  μ                  j                                                      ⁢                          D              ⁡                              (                                  µ                  ,                  y                                )                                              +                                    α              j              n                        ×                          β              j              n                        ×                                          ∂                                                                                              ∂                                  μ                  j                                                      ⁢                                          R                n                            ⁡                              (                µ                )                                                                        Formula        ⁢                                  ⁢                  (          4          )                    
Here αn, βn, and Rn are not direct estimation targets but are auxiliary variables for obtaining a reconstructed image, and are variables for deciding update conditions for obtaining the reconstructed image, which depend on the number of iterations or the pixel position. Such variables are called state variables in this specification. Such state variables may also be called update condition variables in this specification. Here, Rn is a function in a strict sense, but in this specification, those including this will be called state variables.
FIG. 5 is a flowchart showing processing for reconstructing an image by an iterative image reconstruction method using such state variables.
When reconstructing an image by an iterative image reconstruction method, initially, preprocessing is performed. In this preprocessing, the reconstructed image μ0 and the state variables α0, β0, and R0 are initialized (Step S11). At this time, as the initial image, for example, a blank image (image in which all pixel values are zero) may be used.
Next, iterative processing is performed. This iterative processing is a processing step that iteratively performs the calculation represented by the aforementioned Formula (4). At this time, firstly, forward projection processing is performed (Step S12). Next, back projection processing is performed (Step S13). Then, after performing image update processing (Step S14), update processing of state variables is performed (Step S15). At this time, the second term in Formula (4) is calculated by the forward projection processing and the backprojection processing, the third term is calculated by the image update processing, and addition calculations of the entire right side are performed. Then, αn, βn, and Rn are updated. This iterative processing is repeated until the convergence condition is satisfied (Step S16).
Thereafter, post-processing is performed. At the time of this post-processing, the reconstructed image after convergence condition is satisfied is output (Step S17).
The following documents are incorporated by reference in their entirety:    Patent Document 1: Japanese Patent Application Publication No. 2011-156302    Patent Document 2: U.S. Pat. No. 8,958,660    Non-Patent Document 1: C. Lemmens: Suppression of Metal Artifacts in CT Reconstruction Procedure That Combines MAP and Projection Completion IEEE Transactions on Medical Imaging, Volume: 28 Issue: 2 (2009)    Non-Patent Document 2: E. Tanaka: Subset-dependent relaxation in brock-iterative algorithms for image reconstruction in emission tomography, Phys Med Biol 48: 1405-1422 (2003)